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Abstract 

Using molecular dynamic simulations we study a system of particles interacting through a con- 
tinuous core-softened potentials consisting of a hard core, a shoulder at closest distances and an 
attractive well at further distance. We obtain the pressure-temperature phase diagram of of this 
system for various depths of the tunable attractive well. Since this is a two length scales potential, 
density, diffusion and structural anomalies are expected. We show that the effect of increasing 
the attractive interaction between the molecules is to shrink the region in pressure in which the 
density and the diffusion anomalies are present. If the attractive forces are too strong, particle will 
be predominantly in one of the two length scales and no density of diffusion anomaly is observed. 
The structural anomalous region is present for all the cases. 

PACS numbers: 



* Electronic address: 'jonathas@if . ufrgs . brj 
^^Electronic address: _esalcedo@isc._ufsabJ 
Electronic address: oliveir a@iceb . uf op . br | 



I. INTRODUCTION 



The phase behavior of single component systems as particles interacting via the so-called 
core-softened (CS) potentials are receiving a lot of attention recently. These potentials 
exhibit a repulsive core with a softening region with a shoulder or a ramp [JQ. These 
models were motivated by the aim of construct a simple two-body isotropic potential capable 
of describing the complicated features of systems interacting via anisotropic potentials. This 
approach generates models analytically 13l-ll6l| and computationally l|-[9[ tractable still 



capable to retain the qualitative features of the real complex systems. 

The physical motivation behind these studies is the recently acknowledged possibility 
;hat some single component systems display coexistence between two different liquid phases 



17Hl9l|. a low density liquid phase (LDL) and a high density liquid phase (HDL), ending 



at a LDL-HDL critical point. This opened the discussion about the relation between the 
presence of two liquid phases, the existence of thermodynamic anomalies in liquids and the 
form of the potential. The case of water is probably the most intensively studied. A liquid 
where the specific volume at ambient pressure starts to increase when cooled below T ^ 4°C 
, 2l|. Besides, in a certain range of pressures, water also exhibits an anomalous increase 



2G 



of compressibility an d sp ecific 



for Te, 



leat upon cooling from experiments 22|, |23|. Experiments 



M Ga, Bi, |25|L S, [26 



28l | and simulations for silica, 29|-|31| 



and GeisTess, 

silicon [32| and BeF2, ^9|| show that these materials present also density anomaly. 

Besides the anomalies discussed above, water has dynamic anomalies as well. Experiments 
show that the diffusion constant, D, increases on compression at low temperature, T, up to 
a maximum D^i^x{T) at P = PDmax(^)- The behavior of normal liquids, with D decreasing 
on cornpression, is restored in water only at high P, e.g. for P > Pomax ~ 1-1 kbar at 



10°C 



21 



22| . Computational simulations for the Simple Point Charge/Extended (SPC/E) 



water model 



33| recover the experimental results and show that the anomalous behavior 



of D extends to the metastable liquid phase of water at negative pressures - a region that 
is difficult to access for experiments 34-3^. In this region the diffusivity D decreases for 
decreasing p until it reaches a minimum value D^i^^iT) at some pressure PDmm(T), and 
the normal behavior, with D increasing for decreasing p, is reestablished only for P < 



Besides water, silica 



31 



and silicon 



40| also exhibit a diffusion 



anomalous region. 



2 



Acknowledgin K th at CS potentials might engender density and diffusion anomahes, de 
Ohveira et al. 0j41-45| proposed a simple CS model. It has a repulsive core that exhibits 
a region of softening where the slope changes drastically. This model exhibits density, 



diffusion and structural anomalies like the anomalies present in experiments 
simulations 



21 



34 



22| and 



36| for water. This simple system has no attraction between the particles 



and, therefore, no liquid-gas or liquid-liquid critical points are present. Realistic models 
should have attractive interactions since most molecules attract each other either due to van 
der Waals interactions or to more sophisticated electrostatic forces. 

Which effect in the pressure-temperature phase diagram one might expect from the addi- 
tion of a larger attractive part in the potential? For one length scale potentials, the increase 
of the attractive well leads to an increase in the temperature of the liquid-gas critical point. 
In the case of the continuous two length scale potential the same behavior might be expected 
for the liquid-gas critical point but it is not clear which effect the depth of the well has in 
the location in the pressure-temperature phase diagram of the liquid-liquid critical point. 
Moreover, it is also not clear which effect the attraction has in the location in the pressure- 
temperature phase diagram of the density, diffusion and structural anomalous regions. 

In this paper we address these two questions by studying the pressure-temperature phase 
diagram of a potential with a repulsive core followed by a tunable attractive well. We check 
if the introduction of the attraction between particles affects the liquid-liquid critical point 
and the density, diffusion and structural anomalies. 

The remaining of this paper goes as follows. In Sec. Ullthe model is introduced and the 
methods are presented. Details of simulations are given Sec. IIIII In Sec. [IV]the results are 
discussed and, finally, the conclusion are made in Sec. |Vl 

II. THE MODEL 



The model consists of a system of particles of diameter a, inside a cubic box with 
volume V , resulting in a number density p = N/V . The interacting effective potential 
between particles is given by 



U*{r) 
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+ aexp 



r — Tq 



a 



+ feexp 
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r — n 
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where U*{r) = U{r)/e. The first term of Eq. ([T]) is a Lennard- Jones potential of well depth 
e. The second and third terms are Gaussians centered on radius r = ro and r = ri, with 
heights a and b, and widths c and d respectively. This potential can represent a whole family 
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47| 



of two length scales intermolecular interactions, from a deep double wells potential 
to a repulsive shoulder [5i], depending on the choice of the values of the parameters. 

For 6 = the attractive part vanishes and the potential becomes purely repulsive. This 
case was previously studied for determining the pressure-temperature phase diagram as well 
as the regions where water-like anomalies occur 

How the addition of an attractive part in the potential affects the overall pressure- 
temperature phase diagram? In order to answer to this question we obtain the pressure 
temperature phase diagram of the potentials illustrated in Fig. [T] where the attractive part 
is increased systematically without changing the core-softened part of the potential. This is 
done by setting the potential given by Eq. ([1]) with fixed parameters: a = 5.0, ro/o" = 0.7, 
c = 1.0, ri/a = 3.0, d = 0.5 for the five cases studied in this work. The parameter b for 
each case is shown in Table [H 

TABLE I: Parameter b in the potential Eq. ([T|) for each case studied in this work. The other 
parameters are a = 5.0, ro/o" = 0.7, c = 1.0, ri/a = 3.0, and d = 0.5 for the five cases. 
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The potential shown in Fig. [T] has two length scales within a repulsive shoulder followed 
by a attractive well. The addition of an attractive part to the ramp-like format gives rise 
to a liquid-liquid first order phase transition and to a first order liquid-gas phase transition 
ending at critical points. The liquid-liquid phase transition is located in the vicinity of the 
anomalous region. 
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FIG. 1: Interaction potential Eq. ([T]) with parameters a = 5.0, r^/a = 0.7, c = 1.0, ri/a = 3.0 and 
d = 0.5 for all cases, h is shown in Table H] for each case. 

III. DETAILS OF SIMULATIONS 




m 



For the case in which 6 = the results shown in this paper were adapted from Refs. 



|41| |. For the other cases {h ^ 0) the details of simulations are as follows. 
The quantities of interest were obtained by A^V^T-constant molecular dynamics using 



the LAMMPS package 
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= 1372 particles were used into a cubic box with periodic 
boundary conditions in all directions. The interaction through particles, Eq. ([T]), had a cutoff 
of 4.5(7 and the Nose-Hoover heat-bath was used in order to keep fixed the temperature. 

All simulations were initialized in a liquid phase previously equilibrated over 5 x 10^ steps 
at T* = 0.6. The time step used was 0.001 in reduced units and the runs were carried out for 
a total of 3 X 10^ steps, dumping instantaneous configurations for every 2000 steps, giving 
then a total of 1500 independent configurations. The first 200 configuration were discarded 
for equilibration purposes, thus 1300 configurations were used for sampling averages. 
Temperature, pressure, density and diffusion are shown in dimensionless units. 



e 

3 



p = pa 



5 



D-^£^^. (2) 
a 

The pressure of the system is calculated by means of the the virial theorem, 



P = pA:BT + ^^^f(r,,)r,,j, (3) 



where Vij is the vector that it connects particle i with particle j, f(r) =-VU(r). The symbol 
(...) indicates ensemble average. 

The mobility of particles is evaluated by the mean square displacement, given by 

(Ar(r)2) = ([r(ro + r)-r(ro)]^) . (4) 

The diffusion coefficient is then obtained from the expression above by taking the infinite 
time limit, namely 

D=lim^^. (5) 

For normal fluids the diffusion at constant temperature grows with decreasing density. 
Actually in most cases it is expected that it would follows the Stokes- Einstein relation, i.e., 
DocT. 

The structure of the system studied by using the translational order parameter, defined 

as 



31 



36 



49| 



t= f^^m)-Md^, (6) 







where ^ = rp^/^ is the inter-particle distance divided by the average separation between pairs 
of particles p~^^^. g is the distribution function of pairs, is the distance cutoff, where 
we use half of the length of the simulation box, r^, multiplied by p^/^. Another alternative to 
Tc would be the first or the second peak in the g (r). Our choice is preferable, first, because it 



is the maximum distance allowed for the calculation of g (r) 50| giving us a better approach 
allowed for t. Second, the peaks of g (r) change place according to density and temperature 
of the system. Thus additional work would be necessary to find such positions. 

For the ideal gas, g = 1 thus t = 0. As the system becomes more structured a long range 
order {g ^ 1) appears and t assumes large values. The translational order parameter has its 
maximum value in the crystal phase . Therefore, t gives a measurement of how close is the 
fluid close to the crystallization. For a fixed temperature normal fluids present a monotonic 
t{p) curve, increasing with density. 



IV. RESULTS 



Phase Diagram 

Fig. |5] illustrates the pressure-temperature phase diagram for the cases A-E obtained 
through simulations using the potential shown in Fig. [1] As the attractive well becomes 
deeper, the liquid-gas critical point appears and goes to higher temperatures what can be 
easily understood as follows. At low densities the liquid-gas transition is observed by cluster 
expansion namely 

- l-2np [ f{rydr-^^^ [ [ [ f {r)f{r')f{\r - r'\) sin er\''^drdr'd9 (7) 



p ' J ' 3 

where /(r) = e~^^^'^'^ — 1. The critical point is located at 

- 

dp 

The low density behavior obtained using the cluster expansion is illustrated in 
Fig. [3l For T* = 0.60 Fig. |3] shows the pressure-density phase diagram for b = 
0.0, —0.25, —0.50, —0.75, —1.00 using the second and the third virial. For b = —1.00 
the unstable region of the pressure-density phase diagram is large and the system at 
this temperature is deep in the liquid-gas coexistence region of the pressure-temperature 
phase diagram. For b = —0.75 the unstable region is present but is rather small. For 
b = 0.0, —0.25, and —0.50 no unstable region in the pressure-density phase diagram is ob- 
served indicating that the system is above the liquid-gas transition and that T = 0.60 
is larger than the critical point temperature. The comparison between the cases with 
b = 0.0, —0.25, and —0.50 suggests that since the slope of the pressure-density phase di- 
agram increases as b increases, the liquid-gas critical temperature decreases as b increases, 
T*{b = -0.25) < T*{b = -0.50) < T*{b = -0.75) < T*{b = -1.00). Consequently the 



attractive part favors the liquid p 
served in discontinuous potentials 



lase to exists for higher temperatures what is also ob- 



51 



52[ |. Fig. m obtained from the simulations illustrated 
in Fig. ?? summarizes the effect of the attractive part in the location of the critical points 
in the pressure-temperature diagram. 

At high densities where the liquid-liquid phase transition is present the cluster expansion 
with second and third virial is not appropriated. Simulations show that as b decreases 
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FIG. 2: Pressure-temperature phase diagram for the five cases studied in this work. The gray hues 
are the isochores. (a) Case A (6 = 0): p* = 0.04, 0.06, 0.07, 0.08, 0.09, 0.10, 0.107, 0.11, 0.115, 
0.120, 0.125, 0.130, 0.134, 0.140, 0.144, 0.148, 0.154, 0.158, 0.160, 0.168, 0.174, 0.180, 0.188, 0.194, 
and 0.20 from bottom to top. (b) Case B (6 = -0.25): p* = 0.01, 0.015, . . ., and 0.165 from bottom 
to top. (c) Case C (6 = —0.50): same as panel (b), (d) Case D (6 = —0.75): same as panel (b). 
(e) Case E (6 = -1.00): p* = 0.02, 0.025, . . ., and 0.2 from bottom to top. The solid, bold line 
is the TMD line, the dashed line mark the maxima and minima in the diffusion and the dotted 
line bounds the region of structural anomaly. The filled and open circles are the liquid-liquid and 
liquid-gas critical points respectively. 
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FIG. 3: Pressure versus density for T* = 0.60 for the cases b = 0.0, —0.25, —0.5, —0.75, —1.00 from 
top to bottom. 

the pressure needed to form the high density Uquid phase, decreases. This reflects that 
the attractive part favors the high density Uquid phase over the low density liquid phase. 
The attraction leads in this case to a more compact liquid phase what is also observed in 



discontinuous potentials |5ll. |52|. 



Density anomaly 

The density anomalous region in the pressure temperature phase-diagram can be found 
as follows. From the Maxwell relation, 

the condition for density anomaly at fixed pressure, i.e., a maximum in p(T) curve given by 
{dp/dT)p = 0, is equivalent to the condition (dP/dT)^ = 0, corresponding to a minimum in 
the -P(T) function. While the former is suitable for A^PT-constant experiments/simulations 
the latter is more convenient for our A^V^T-ensemble study, thus adopted in this work. 

In this sense, the minima at the isochores mean that the system has density anomaly. 
These extrema points in the pressure-temperature phase diagram are named as temperature 
of maximum density (TMD) points, which connected form the TMD line. Fig. [2] shows 
the TMD line as a solid, bold line in panels (a)-(d), corresponding to the cases in which 
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FIG. 4: Location of the critical points for the cases B-E considered in this work. The Case A does 
not present any fluid-fluid critical point whereas the Case B has a liquid-gas but no liquid-liquid 
critical point. The symbols have the same meaning as in Fig. [21 i.e., filled and open circles mark 
the liquid-liquid and liquid-gas critical points respectively. The arrows indicate the direction of 
increasing the attractive interaction. 



b = 0.0, —0.25, —0.50, and —0.75 respectively. As the attractive well becomes deeper, the 
region in the pressure temperature phase-diagram occupied by the density anomalous region 
shrinks and moves to lower pressures and higher temperatures until to the limiting case 
[b = —1.00) in which no density anomaly is present [note that there are no local minima in 
the P(T) curves of Fig. [SJ^e)]. 

The link between the depth of the attractive region and the presence or not of the 
TMD goes as follows. The TMD is related to the presence of large regions in the system 



in which particles are in two pre: 



second scale in our potential 12 



'erential distances represented by the first scale and the 
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55| . While for normal liquids as the temperature is 



increased the percentage of particles at closest scales decreases [see the case (e) in Fig. \5\, 
for anomalous liquids [see cases (a), (6), (c) and (d) in the Fig. [5] there is a region in the 
pressure-temperature phase diagram where as the temperature is increased the percentage 
of particles at the closest distance increases while the percentage of particles in the second 
scale decreases. This increasing in the percentage is only possible if particles move from the 
second to the first scale. 

In Fig. El^e), the decrease of particles in the first scale leads to a decrease of density with 
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FIG. 5: Radial distribution function versus distance for (a) the Case A (6 = 0.0) with p* = 0.14 
and T* = 0.25,0.35,0.45,0.55,1.0,2.0,3.0, and 4.0; (b) Case B (6 = -0.25) with p* = 0.085 
and T* = 0.32,0.36,0.44,0.56,0.68, and 0.80; (c) Case C (6 = -0.50) with p* = 0.06 and 
T* = 0.44,0.48,0.60,0.72,1.0,1.6,2.0, and 3.5; (d) Case D (6 = -0.75) with p* = 0.06 and 
T* = 0.40, 0.44, 0.48, 0.56, 0.60, 0.68, 0.72, 0.80, and 1.0; and (e) Case E (5 = -1.00) with p* = 0.07 
and T* = 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, and 0.80. The arrows indicate the direction of in- 
creasing temperature. 
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an increasing temperature: behavior expected for normal liquids. In the Fig. [5]^a)-(d), the 
increase of particles in the first scale leads to an increase of density with temperature what 
characterizes the anomalous region. The anomaly is, therefore, related to the increase of 
the probability of particles to be in the first scale when the temperature is increased while 
the percentage of particles in the second scale decreases. As the potential becomes highly 
attractive this "mobility" between scales disappears, i.e., the high density liquid becomes 
dominant and no anomalous region is observed. 

Diffusion anomaly 

The mobility of any liquid is given by the diffusion constant. Figure |6] shows the behavior 
of the dimensionless translational diffusion coefficient, D*, as function of the dimensionless 
density, p*, at constant temperature for h = 0.0, —0.25, —0.50, —0.75, and b = —1.00. The 
solid lines are polynomial fits to the data obtained through simulation (dots in Fig. [6]). For 
normal liquids, the diffusion coefficient at constant temperature decreases with density. For 
the cases A-D [shown in Fig. [6](a)-(d)] D* anomalously increases with density in a certain 
range of pressures and temperatures instead. From Figure [6] we see that for very small and 
very high densities D* decreases with increasing density as expected for a normal liquid. 
For intermediate values of density, pomax > P > PDmin, D* increases with increasing density 
what leads to local maxima at pomax and a local minima at pDmin- These local extrema 
in the diffusion versus density plots bound the region inside which the diffusion behaves 
anomalously (dashed lines in Fig. [6]). This region is mapped into the pressure-temperature 
diagram illustrated in Fig. [2] as dashed lines in (a)-(d). As the attractive well becomes 
deeper, the diffusion anomalous region in the pressure-temperature phase diagram shrinks 
and it goes to lower pressures. In the case in which h = —1.00, shown in Figure El^e), the 
diffusion constant behaves as in a normal liquid. This result again is consistent with the 
idea that a deeper attractive term favors the high density liquid phase. 

Structural anomaly 

Besides the density and the diffusion anomalies an structural anomalous region might be 
present. Figure [7] shows the translational order parameter as a function of density for fixed 



12 



0.4 


.^^^ (a) _ 










0.3 












^ 








1 \ 




0.2 




1 1 




0.1 


\ "^-^^—"""^ 1 ^^--^^ 


1,1 




0.06 0.08 0.1 0.12 




0.04 0.06 0.08 0.1 




0.04 0.06 0.08 0.1 

* 

P 



0.04 0.08 0.12 



FIG. 6: The diffusion coefficient against density for the (a) Case A, with isotherms 0.2, 0.23, 0.262, 
0.3, 0.35, 0.4, and 0.45 from bottom to top. (b) Case B with isotherms 0.16, 0.20, . . ., and 0.56, 
(c) Case C, whose temperatures shown are 0.36, 0.40, . . ., and 0.68, (d) case D, with temperatures 
0.48, 0.52, and 0.80, and (e) Case E with isotherms 0.70, 0.75, 1.0, 1.10, and 1.70. 
The dashed hues mark the local maxima/minima in the D{p) curves. For the region enclosed by 
these lines particles move faster under compression. The dashed lines in this figure have the same 
meaning as those ones in Fig. [2j 
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temperatures for the potential we are studying for b = 0.0, —0.25, —0.50, —0.75, and —1.00. 
The dots represent the simulation data and the solid lines are polynomial fit to the data. 

The non-monotonic behaviour of these curves indicate that there is a region in which t 
decreases with density. This means that the system becomes less structured for increasing 
density. Dotted lines determine the local maxima and minima of t, bounding the structural 
anomalous region. This region was mapped into the pressure-temperature phase diagram 
(dotted lines), as can be seen in Figure |2l The comparison between the behavior for different 
b values indicates that as the attractive well becomes deeper the structural anomalous region 
in the pressure-temperature phase diagram shrinks and moves to lower pressures and it is 
still present even in the deepest case, b = —1.00. According to these results we believe that 
for b < —1.00, i.e., cases in which the attractive part is more intense than one showed in 
Case E, the structural anomalous region will also vanish. This result again is consistent 
with the idea that a deeper attractive term favors the high density liquid phase. 

Figure M gives an overview of the density, diffusion, and structural anomaly locations in 
the pressure-temperature phase diagram. 

V. CONCLUSIONS 

In this paper we have explored the effect of the addition of an attractive part in a two 
length scales potential. Particularly we analyze if the depth of the attractive part changes 
the position (and the presence of not) in the pressure-temperature phase diagram of the 
two liquid-gas and liquid-liquid critical points and of the density, diffusion and structural 
anomalous regions. 

For sufficiently intense attraction between particles both the liquid-liquid and the liquid- 
gas critical points are present. These two critical points are observed even for a very at- 
tractive potential. For a small attractive interaction, only the liquid-gas critical point was 
found what indicates that for the coexistence of two liquid phases the attractive well have 
to be deeper than a certain threshold. 

Since the attraction favors the liquid phase (particularly the high density liquid phase), as 
the b decreases the liquid-gas critical point moves to higher temperatures (shown in Fig. H]) 
and the liquid-liquid critical point to lower pressures. 

The density, diffusion and structural anomalous regions are present even in the absence of 
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FIG. 7: Translational order parameter against density for (a) Case A, where each hne correspond 
to an isotherm. The isotherms are: 0.25, 0.30, . . ., 0.55, 0.7, 1.0, 1.5, 2.0, and 2.5 from top to 
bottom, (b) Case B, with isotherms 0.20, 0.28, . . ., 0.68, 0.80, 1.0, 1.2, 1.6, 2.0, and 2.5 from top 
to bottom, (c) Case C whose temperatures are 0.36, 0.40, . . ., 0.80, 1.0, 1.2, 1.6, 2.0, 2.5, and 3.0 
from top to bottom, (d) case D with T* = 0.52, 0.56, . . . , 0.80, 1.0, 1.2, 1.6, 2.0, 2.5, 3.0, and 3.5 
from top to bottom. Finally, (e) case E with T* = 0.70, 0.75, . . ., 1.0, 1.10, . . .,1.70, 2.0, 2.5, and 
3.0 from top to bottom. The dotted lines bound the region of structural anomalies, i.e., the region 
where the parameter t decreases upon increasing density. 
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FIG. 8: (a) TMD line for the cases considered in this paper. Note that there is no density anomaly 
in the Case E (see Fig. [2]). (b) The diffusion anomaly region for the Case A-D. No diffusion 
anomaly was found for the Case E (see fig. (H). The shadowed regions correspond to the region 
between the dashed lines in Fig. [2l In (c) is shown the structural anomalous region for Case A-E. 
Here, the shadowed region corresponds to the region between the dotted lines in Fig. [5J See the 
text for discussion. 

attraction. As b decreases, the high density liquid structure is favored and so the anomalous 
regions in the pressure-temperature phase diagram ( shown in Fig. [8]) shrinks, moves to 
lower pressures and disappears for very attractive potentials. 

In resume density and diffusion anomalous regions are present in two length scales po- 
tential if the attractive interaction is not too strong. 
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